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Abstract 

The method of generahzed modehng has been apphed successfully in many different 
contexts, particularly in ecology and systems biology. It can be used to analyze the 
stability and bifurcations of steady-state solutions. Although many dynamical systems 
in mathematical biology exhibit steady-state behaviour one also wants to understand 
nonlocal dynamics beyond equilibrium points. In this paper we analyze predator- 
prey dynamical systems and extend the method of generalized models to periodic 
solutions. First, we adapt the equilibrium generalized modeling approach and compute 
the unique Floquet multiplier of the periodic solution which depends upon so-called 
generalized elasticity and scale functions. We prove that these functions also have to 
satisfy a flow on parameter (or moduli) space. Then we use Fourier analysis to provide 
computable conditions for stability and the moduli space flow. The final stability 
analysis reduces to two discrete convolutions which can be interpreted to understand 
when the predator-prey system is stable and what factors enhance or prohibit stable 
oscillatory behaviour. Finally, we provide a sampling algorithm for parameter space 
based on nonlinear optimization and the Fast Fourier Transform which enables us to 
gain a statistical understanding of the stability properties of periodic predator-prey 
dynamics. 

Keywords: Generalized models, periodic orbits, predator-prey system, Floquet theory, 
moduli space flow, Fourier series, discrete convolution, parameter sampling, optimization, 
correlation. 



1 Introduction 

Predator-prey systems have been a cornerstone in mathematical biology for many decades 
[4]. Standard textbooks on dynamical systems, differential equations and ecology provide 
a plethora of models that aim at capturing the interaction between a predator population 
Y and a prey population X. Examples for modeling the situation by ordinary differential 
equations (ODEs) are |6l |5] 
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where Pi, ki are parameters and (LV) are the Lotka-Volterra equations and (RM) is the 
Rosenzweig-MacArthur model |32J. These two models are the most common examples of a 
large class of different models of the form 



where a > is a parameter describing biomass conversion efficiency and the functions S, G 
and M represent prey growth, predation, and predator mortality, respectively. Because the 
parameter a can always be removed by scaling the variable Y and re-labelling the functions 
we will always assume a = 1 from now on. 

Generalized models [201 ES] directly work with the formulation ([1]) without specifying 
functional forms for G and M. Previous works on generalized models [HI [151 [I2l [52] 
focused on analyzing the dynamics close to stationary states. Beyond the structure of the 
equations ([T|) this analysis requires only the assumption that steady states exist in the class 
of models under consideration. The central idea of generalized modeling is to parametrize 
all possible Jacobians that can be encountered in steady states in the class of systems under 
consideration. Using a specific renormalization procedure, one can define parameters that are 
easily interpretable (and often also directly measurable) in the context of an application. Ap- 
plications of generahzed models to ecology can be found in [T^l [TSj [TTj I2U1 [^ 1^ liTf 12^ liSl I5U] . 
Let us emphasize that we do not claim that stability results from generalized models have 
not been observed before in some specific models; in fact, the literature on stability of planar- 
predator prey systems is very large. For instance, questions of local and global stability have 
been investigated in various predator-prey systems [131 [211 [IS [IHl [36] . 

In the present paper we go beyond the previous analysis and study nonstationary dy- 
namics in the context of generalized modeling. We extend the theory of generalized models 
to arbitrary periodic solutions in the context of the predator-prey system ([1]). We show that 
this mathematical extension of generalized models yields several new phenomena in com- 
parison to generalized models for steady states. For example, we can define time-periodic 
generalized parameters of the predator-prey model and we prove that these functions obey 
a system of ODEs (a flow on moduli space). Using Floquet theory [7] and Fourier analysis 
[28] we derive analytical conditions for the solvability of the moduli space flow and obtain 
an analytical stability formula. In this context, a main result is that the stability formula 
for periodic solutions only depends on two constants that can be calculated via a discrete 
convolution. Using this formula we can identify parameters and conditions that enhance 
the stability of predator-prey cycles. Furthermore, we develop an algorithmic approach to 
sample the function space of parameters by solving an auxiliary optimization problem, which 
will be instrumental for future applications to larger systems. 

The paper is structured as follows: In Section [21 we recall the necessary tools from steady- 
state generalized models, Floquet theory, and Fourier analysis. In Section [21 we calculate 
the generalized vector field for non-equilibrium solutions. In Section [H we derive the flow 
on moduli space. In Section [5l we compute the generalized scale and elasticity functions for 
several specific functional forms to gain a better understanding how generalized and specific 
models link up. In Section [HI we use tools from Fourier analysis to derive algebraic conditions 
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from tlie moduli space flow. In Section [71 we provide an analytical stability analysis of 
periodic orbits in generalized predator-prey models. Using this result we can identify which 
situations increase or decrease stability and interpret the results in an ecological context. 
In Section [HI we develop a sampling technique for generalized scale and elasticity functions 
that is based on solving an auxiliary optimization problem and the Fast Fourier Transform. 
We also use this sampling-based approach to improve our understanding of stabilizing and 
destabilizing factors of the predator-prey system. In Section [HI we conclude with a brief 
summary and outline the large range of applications and theoretical challenges that can be 
found in non-equilibrium generalized models. 



2 Background 

In this section we introduce essential tools and techniques that will be used throughout this 
work. Further, we use the opportunity to fix the notation. Below, we denote a general 
ordinary differential equation ODE by 

^ = Z' = F{Z), forZeM^ (2) 

and assume always that F is sufficiently smooth. In the following we are going to recall the 
necessary tools from steady state generalized models, Floquet theory and Fourier analysis. 



2.1 Generalized Models 

Let us start by reviewing generalized modeling [20] for ODEs with equilibrium points. A 
detailed mathematical approach to generalized models can be found in [33]. For the present 
discussion we restrict ourselves to review generalized models in the context of a planar 
predator-prey system [21]. Such systems describe the interaction of a population of prey X 
and a population of predators Y. The prey population grows at rate S{X), predation occurs 
at rate G{X,Y) and natural mortality of the predator at rate M{Y). Denoting the prey 
density as X and predator density as Y we capture the dynamics by 

X' = S{X)-G{X,Y), 
Y' = G{X,Y) - M{Y), 

where S, M e C"^(M+,M+) and G G C"'(M+ x M+,M+) are sufficiently smooth functions. 
Generalized modeling assumes that (|3]) admits an equilibrium point (X, Y) = {X*, Y*) G 
M+ X M+. 



We normalize the equilibrium defining new coordinates 

^ - ^ and y := — . (4) 

This transformation moves the equilibrium to {x,y) = (1, 1). The next step is to normalize 
the rate functions 

s(x)-^^ ,(.,)- ^(^^"'^^'^) m(y)-^^^ (5) 
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A direct substitution of dl])-© into ([3]) gives 



/ s(x*) I N G(x*,y*) / X 



where the prefactors of the form S{X*) / X* ^ G{X* ^ Y*)/ X* , etc. represent normahzed fluxes 
in the steady state and are also called scale parameters 

_S{X*) _ G{X\Y*) _ G{X*,Y*) _ M{Y*) 

Ps ■- , Pi-- J^^ , P2 ■- , Pm ■- ■ [') 

Since {x,y) = (1,1) is an equilibrium point we know that the following holds: 

= ^^^(l,l)-^m(l) = 



Therefore (ED can be re-written as 



x' = /3i{s{x) - g{x,y)), 
y' = I32i9ix,y) - m{y)). 

The Jacobian at the equilibrium {x,y) = (1, 1) is then given by 



I329x P2[9y-'my\ 
where dx, dy denote partial derivatives and we refer to the constants 

Sx = dxis{x))\x=i, 9x = dxi9ix,y))\(x,y)={i,i), 

9y = dy{9{x,y))\^x,y)={i,i), rriy = dy{m{y))\y=i, 



(9) 



jn I) = I (^1 9x[six) - g{x,y)]\^x,y)={i,i) -f^i dy[g{x,y)]\(^x,y)=(i,i) ^ 
' 1^2 dx[9ix,y)]\(^x,y)=ii,i) 1^2 dy[g{x,y) -m{y)]\^^^y)=i^i ' 



:iii 



(12) 



as elasticities. The scale parameters and elasticities are also referred to as generalized pa- 
rameters. 



In the following we will use the insight that every power law function corresponds to an 
elasticity that is identical to the exponent of the power law. For example, if we assume that 
the mortality M{Y) is a linear function M{Y) = KY then we find 

( M{Y*y) \ fKY*y\ 

Hence we can relate the growth properties of the unspecified functions forms to the elastici- 
ties. 

The stability of the equilibrium {x,y) = (1,1) can be inferred from the eigenvalues of 
J(l, 1) and hence only depends on the generalized parameters. This admits a bifurcation 
analysis of all steady state models of the form in generalized parameter space. Despite 
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the large class of models that one treats simultaneously it is often easy to interpret scale 
parameters and elasticities in applications [20] • Thereby a generalized model enables us 
to draw conclusions about a whole class of differential equations, for further examples see 

[221 SHI ig mug. 

We note that generalized modeling can also be applied to equilibria for delay equations 
[271 [33], spatially homogeneous states for partial differential equations [3] and to stochastic 
differential equations [53] . 

2.2 Floquet Theory 

For analyzing the stability of periodic solutions in GM we resort to the framework offered 
by Floquet Theory. Suppose has a period orbit ■jit) = jit + T) with minimal period 
T. Let S denote a suitable (A^ — l)-dimensional transversal section to F and consider the 
associated Poincare map P : S — )■ S. This map has a fixed point C S associated to 
the periodic orbit 7 i.e. -P(X^) = X^. Recall [HI El] that the stability of 7 is determined 
by the — 1 eigenvalues (or characteristic/Floquet multipliers) Ai, . . . , Aat-i of the matrix 
DP{X^). If |Aj| < 1 for all j G {1, . . . , — 1} then the periodic orbit is stable, if there 
exists Xj such that \Xj\ > 1 then the orbit is unstable and eigenvalues with \Xj\ = 1 signal 
bifurcations under parameter variation. We can study the stability of 7 by considering the 
non-autonomous linear variational equation 

v' = DF{j{t))v =: A{t)v (13) 

where A{t) is periodic. An N x N matrix M{t) that satisfies 

M' = A{t)M with M(0) = Id (14) 

is called the fundamental matrix solution of f ll3p . The constant matrix M{T) is called the 
monodromy (or circuit) matrix. It has eigenvalues 

I5 -^17 -^27 • • • ) Xn-1 

where the trivial eigenvalue 1 is associated to the direction tangent to the periodic orbit that 
links the variational equation to the Poincare map P. Furthermore, the Liouville formula 

A1A2 ■ ■ ■ Aiv-i = det M(T) = exp Tr{A{t))d?j (15) 

holds. Floquet's theorem states that there exists a T-periodic coordinate change C{t) and a 
constant matrix R such that 

M{t) = C(t)e*^. 

Since M(0) = Id it follows that C(0) = C{T) = Id and we find that the monodromy matrix 
can be expressed as 

M(T) = e^^. 

An elegant explicit formula for the Floquet multiplier from ( !T5|) is only available for N = 2. 
In general the computation of Floquet multiplier thus requires numerical approaches, which 
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typically start with computing the periodic solution with a suitable boundary value method 
such as collocation or finite differences |3H E]- The variational equation (HM is solved on 
suitable sub-intervals of the periodic orbit discretization as an initial value problem to ob- 
tain M(T). The eigenvalues of M(T) are then obtained yielding the Floquet multipliers. 
Although, in certain circumstances, such as large multipliers, the computation can be nu- 
merically problematic [101 EZ] • 

Let us point out that Floquet theory has not been widely applied in the context of 
ecology [30j although it is a standard tool in the mathematical theory of dynamical systems 
[8j . Klausmeier [30] suggests that "Floquet theory [is] a useful tool for studying the effects of 
temporal variability on ecological system". In the context of our approach, Floquet theory 
is not only a tool for a particular model but we will also show that it nicely extends to 
generalized models. 



2.3 Fourier Series 

Since we work with periodic solutions to ODEs and also other time-dependent periodic 
functions we briefly recall basic facts about Fourier series to flx normalization constants and 
notation. Assume that / : M — M is T-periodic so that we can identify the domain of / 
as the circle M/(TZ) = S^. We can formally write the complex Fourier series J-'[f] of / as 
follows: 

nfm= E /»exp(^) (16) 



where the Fourier coefficients f{k) are 



I r . f 2mks 



f{k) = - / /(.)exp( —]ds. 



Observe that f{k) = f{—k), where the overbar denotes complex conjugation. Further, 
/(O) = 1/T Jq f{t)dt is the time average of the periodic function. The convergence question 
J^[f]{t) —7- f{t) is extremely intricate depending on the properties of / [531 [2H]- In the 
following, all functions we are going to approximate by Fourier series will be in C""(5'^,M) 
for some sufficiently large r or even r = oo. In this case, uniform convergence is immediate. 
A very important practical result in this context is to control the Fourier coefficients. 

Theorem 2.1 (see [28]). If f e C"^(5\M) then \ f{k) \ = 0{k-'') as \k\ oo. 

Theorem 12. II is a version of the Riemann-Lebesgue Lemma for smooth functions and can 
provide an extremely rapid decay of the Fourier coefficients. This justifies (for the smooth 
case!) dropping higher-order terms \k\ > n for some rather small suitable k G N. The 
remaining sum is expected to be a good approximation to the original periodic function /. 
We write 
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We remark that it can be convenient to re-write the complex Fourier series (fT6|) as a real 
Fourier series 

oo 
k=l 

where the real Fourier coefficients relate to the complex ones by 

f{k) = i(afc - ibk) and f{-k) = ^(a^ + ibk) 

for k &Nq. Another important tool in Fourier analysis we will need are convolutions. Recall 
that the discrete convolution of two periodic functions / and g is defined as 

oo 

{f*g){n)= J2 Hmn-k). 

k=—oo 

Obviously the convolution operator is associate, commutative and distributive. 



Ofc cos 



+ hk sin 



(^) 



3 Non-Equilibrium Planar Predator-Prey Systems 

We return to the planar predator-prey system ([3]) from Section (12. ip given by 

X' = S{X)-G{X,Y), 

Y' = G{X,Y) - M{Y). ^ ' 

Denote the vector field of f|T7|) by F{X,Y). The vector field is only considered on the 
ffist (positive) quadrant F : R"*" x M+ — > as predator-prey densities are assumed to be 
non- negative. 

We want to analyze the class of vector fields f[T7l) under the assumption that it admits 
a non-equilibrium orbit that is bounded as |t| — )■ oo. From an ecological point of view the 



most interesting case are limit cycles, so-called predator-prey cycles. We assume that ([E 
has a periodic orbit 7(t) = (7i(t), 72(i)) with period T. The definition of the model implies 
that 7i > for i G {1,2} and all t. 

In the following, we are going to slightly extend the notation employed already in Section 
(12. ip by re-using names for variables and generalized parameters. As in the case of equilibria 
one can consider a normalizing coordinate change 

and , I 

7l 72 

which maps the periodic orbit to the point {x,y) = (1,1) =: 1. The ODEs f[T7P and the 
product rule imply 

X' = x'ji + = x'7i + xFi(7) = S{x'ji) - G{x'ji, 2/72), 
Y' = + 2/72 = y'l2 + 2/^^2(7) = G'(x7i, 2/72) - M{y^2)- 
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Therefore the new equations can be written as 



= J-(5(a;7i)-G(x7i,l/72)-a;(^(7i)-G(7))), 
y' = -^(G(x7i,2/72)-M(y72)-l/i^2(7)) ^ ' 

= I (G(x7i, yi2) - M{y^2) - y(G(7) - M(72))) • 

In analogy to the equihbrium case we introduce normahzed functions 

jW:=^ip^. (19) 

5(7i) G(7i,72) M(72) 



and define the scale parameters 



(20) 



which are now time-dependent T-periodic scale functions. We will often suppress the time- 
dependence in the notation and just write, for instance, (3s instead of f3s(t). Using f|T9|) - fl20|) 
in ( !T8|) we find 

x' = I3s[s{x) - x] - I3i[g{x,y) - x], 
y' = Mg{x,y) -y] - Pmlmiy) -y]. 



For applying Floquet theory we linearize ( 1211) around the limit cycle which yields the matrix 

An-f\ - ( - 1] - PA{d.9){l) - 1] -{dyg){l) 

^ ' ^" V (9^9m P,[{dyg){l) - I] - P^[{dym){l) - I] 



l3s{t)[sS) - 1] - l3i{t)[gS) - 1]) ~gy{t) 

9.{t) P2{t)[gyit) - 1] - Am(t)K(t) - 1]) 



We can re-write A{1; t) in terms of the more familiar elasticities, leading to 

A{l;t) = 

where the four time-dependent T-periodic elasticity functions are 

sAt) ■■= 9.{t) := {d,g){l), gy{t) := {dyg){l), m,(t) := (9,m)(l). 

The periodicity and time- dependence becomes more apparent once we write out the detailed 
definitions, for example 

The previous calculations show that we can introduce replacements for the generalized pa- 
rameters for equilibrium points in the context of periodic orbits. In particular, the scale 
parameters and elasticities become time-dependent and periodic. The term "generalized 
functions" is already used in a different context [51j. Therefore, we refer to elasticity func- 
tions and scale functions directly. To analyze the stability of the periodic solution we use 
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Floquet theory (see Section |2^ . For planar systems the stabihty of the periodic orbit is de- 
termined by computing the only non-trivial Floquet multiplier A. Liouville's formula implies 
that 

A = exp (^J^ TT{A{l;t))dt^ 

= exp(^^ f3s{s, - 1) - f3^{g, - 1) + f32{gy - 1) - f3Umy - l)dty (22) 

We can thus express the Floquet multiplier as a function depending on elasticity and scale 
functions. This is analogous to writing the eigenvalues of the Jacobian as functions of the 
generalized parameters in the equilibrium case. 

4 The Moduli Space Flow 

In analogy to the generalized exploration of local dynamics, the stability of the limit cycle 
can be studied by assuming plausible values for the generalized parameters (here, scale and 
elasticity functions). The value of generalized models lies in their ability to cover the whole 
range of possibilities that are plausible in the system. For an unbiased analysis it is essential 
that we consider only those values of parameters that are consistent with the set up of the 
system. For instance, in case of equilibrium generalized models we must demand that the 
parameter values which we assume do not preclude the existence of an equilibrium solution in 
the class of systems. Likewise, only those scale and elasticity functions should be considered 
which are mutually consistent and thus could arise in at least one example system in the 
class of models under consideration. To understand this problem we briefly go back to the 
equilibrium scenario (see Section [2?T1) . Suppose we just choose a set of generalized parameters 

/3i = P2 = 1^*2, = si, = gl, gy = g*, ruy = m*y, (23) 

where we assume that all parameters are positive. One natural question is if there exist 
specific functions S, G and M that lead to the generalized parameters fl23l) . 

Proposition 4.1. Suppose (1231) are given positive generalized parameters. Then there exist 
functions S, G, M and an equilibrium {X,Y) = {X*,Y*) for so that ([7]) and ([12]) hold 
i.e. there exists a differential equation of the form ([3]) that has the given set of generalized 
parameters. 

Proof Pick M(Y) = piY'^v for some pi e M+ then rriy = dy{M{Y*y) /M{Y*))\y=i = m*y. 
Similarly we pick S{X) = p2X^^ and obtain Sx = s*. Using a slight modification of this 
approach we define G{X,Y) = X^^Y^y and get gy = g* as well as gy = g*. We also must 
have Ps = Pi = Pi and Pm = P2 = P2 which translates into the conditions 

p, p,{x*y-~^ (x*)9--i(y*)9«, 
/32^=^ P2{Y*)"'y-^ ^='^ {x*y'-{Y*yy-\ 

We can always choose pi and p2 to satisfy (C2) and (C4). Then we can use X* and Y* to 
satisfy (CI) and (C3). The result follows. □ 
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Although there are certainly many other ways for constructing functions that are consis- 
tent with a given set of scale and elasticity parameters, already the existence of one such set 
of functions proves that the assumed parameter values could be encountered in the class of 
models under consideration. This observation is of central importance for sampling proce- 
dures, by which high dimensional generalized models are typically analyzed ^61 122] . 

For non-equilibrium systems the situation is different since one has to ask whether a 
whole set of given functions 

13 sit), Pmit), /32(i), s^it), g^{t), Qyit), myit), 

can potentially arise from a system of the form ([1] 



Theorem 4.2. Suppose we are given elasticity functions s^, friy, and gy then the scale 
functions have to satisfy the following set of ODEs 

(24) 



P'm = /3m(/32 - /3m)("^y - 1), 



P'2 = P2m-(3m)gy-W2-Pm) + {Ps-f3l)g.). 

Proof. We start by deriving the equation for f3s. We know that f3s = S'(7i)/7i and direct 
differentiation with respect to time via the quotient and chain rules gives 

7i7i'^'(7i) - '^(7i)7i 



(7i)^ 

7i7i'5'(7i)5(7i) 5(7i)7i 



(7i)^5(7i) (7i)^ 

Noting that = 71 S" (71)/ 5 (71) and using the definition of f3s the equation transforms to 

Sxl'iS{-fi) 5(71)7; 



(71)^ (7 



\2 



= £f2iA _ M = _ 1)21 . (25) 

7i 7i 7i 

Since (71,72) is a trajectory of (fTTI) we must have 7^ = »S'(7i) — G(7i,72)- This implies upon 
substitution into (1251) that 

71 

which is the first equation in ( l24l) . The calculation for is similar. For P[ we find 

7i[G'x(7i, 72)7; + (71, 72)72] -7^^(71,72) 



7iG'x-(7i, 72)71 G'(7i, 72) ^ 71 Cj/ (71, 72) 72 (71, 72) 72 _ 7lG'(7i,72) 
G'(7i,72)(7i)^ (71)^ ^(71, 72)72 (71)^ 

9.Pi- + 9yPi-- Pi- 
ll 72 7i 

Pl{9M - Pi) + 9y{P2 - Pm) - iPs - Pi)) 
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The calculation for fS!^ is similar to the one for I3[ 



□ 



The main conclusion is that the elasticities and scale functions which parametrize the 
ODE fl2T|) satisfy an ODE themselves. Because one often uses the terms "parameters" and 
"moduli" interchangeably, Theorem 14. 2l implies that the time-dependent parameters of gener- 
alized models generate a flow on moduli space. The following remark describes the relevance 
of this viewpoint in some other research areas. 

Remark: The term "moduli space" is perhaps most commonly used in algebraic geom- 
etry which, broadly speaking, is the study of solutions of algebraic equations [211 [25] . The 
solutions form algebraic varieties (e.g. curves). Often suitable parametrized families of al- 
gebraic varieties again have the structure of an algebraic variety, where the latter object is 
the moduli space of parametrized families. The study of the geometry of moduli spaces has 
also been transported into different branches of physics such as quantum field theory In 
dynamical systems theory, a classical moduli space argument is made in the renormaliza- 
tion analysis of parametrized families one- dimensional maps [23] , where the renormalization 
transformation can be viewed as a map generating a dynamical system on moduli space. A 
very similar situation occurs for billiard dynamics where the so-called Teichmiiller flow on 
the space of lattices appears [HI 138] . 

We note that the positive quadrant is an invariant set for (124]) which means that this 
property lifts from the predator-prey family of vector fields to the moduli space. From 
Theorem 14.21 we can immediately infer a condition for the existence of a generalized model 
with given elasticities. 

Corollary 4.3. Suppose Sx, my, and Qy are given T-periodic elasticity functions with 
minimal period T . If flM]) has no T-periodic solutions then there exists no generalized model 
of the form ( fTTl) for the given elasticities. 

Note that the existence of periodic solutions in Corollary 14. 31 is only a necessary condition 
for the existence of a generalized model. We also observe that for the equilibrium case the 
conditions /3i = Ps and (32 = Pm give /3g = = = = consistent with steady state 
generalized modeling. It is also interesting to ask what happens if we do not specify the 
elasticities. 

Taking the idea of deriving a differential equation one step further we consider 




^(7i)[7;5'(7i) + 7i7l5"(7i)] - 7i5'(7i)7i'5'(7i) 





Applying similar substitutions as in the proof of Theorem 14.21 we obtain 




^(7i) [71^^(71) +7171^^^71)] -71^^(71)71^X71) 
l[ , 7;5"(7i) l[ 2 
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Similar calculations can be carried out for g'^, g'y and m'y. These suggest that specifying a 
suitable scaled version of second partial derivatives of S, G and M will provide a system of 
eight ODEs. This procedure could be continued iteratively. It is interesting to note that 
closing a system of ODEs at a given order is a problem that also occurs in the context of 
moment closure for networks [291 HB] and for moment equations of stochastic differential 
equations [15| [TT] . To get a better understanding of the stability of non-equilibrium gen- 
eralized models and the flow on moduli space we proceed to consider a few typical specific 
functions S, G and M that appear in predator-prey models. 



5 Specific Functions 

In this section we calculate the generalized elasticity and scale functions for several well- 
known predator-prey models. All the model parameters ki (for / G N) we are going to use 
below are positive due to modeling considerations. We start with the growth of the prey 
S{X). Typical choices are 

S{X) = kiX (linear growth), 

S{X) = kiX'P (power growth), 

s\x) = kiX-kiX^ (logistic growth), 

S{X) = kiX{k2 -X){X - hi) (growth with strong Allee effect), < A;2 < k^. 



We start by looking at linear growth. We find 



— — ki, Sx — dx ( ^, s 

71 V Shv 



x=l 



For estimating the impact of linear growth on stability we consider the formula ( l22l) and 
view it as a product of exponentials. The term involving fSg and Sx is 

exp(^^ (3,[sx-l]dty (26) 

Therefore, a linear prey growth does not contribute to the non-trivial Floquet multiplier 
because Sx = ^ and exp(Jg Odt) = 1. Different types of polynomial growth with a single 
term can be treated analogously since for S{X) = kiX'^ we find 



Ps = = /i;i7i ) Sx = dx ' 



71 



V ^(7i) 



= p. 

x=l 



where the elasticity function coincides with the result for equilibrium generalized models. 
This allows us to write 

/3,[s.-l] = A;i7r'b-l]. 

Considering ( l26l) we find that increasing p increases the Floquet multiplier and therefore 
has always a destabilizing effect, whereas decreasing p has a stabilizing effect. For logistic 
growth we obtain 



Ps = = ki- /C271, Sx = Ox . 

7i V 7i 



2+ 



x=l 



ki + k2'ji 
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This implies 

/3s[s^ - 1] = ih - A;27i) 
Considering fl26l) again we find 



2 + 



-ki + A;27i 



(27) 



exp 



T 



l]dt 



exp 



exp 



where the integral is positive because k2 > 0. This means that increasing k2 or increasing 
/q 7i(it will promote stability as the Floquet multipher will move closer to 0. For logistic 
growth increasing k2 corresponds to decreasing the carrying capacity ki/k2 of the population. 
This can be interpreted as a manifestation of the paradox of enrichment [13| [T7] which 
captures the observation that increasing the carrying capacity generally has a destabilizing 
effect on attractors observed in ecological models. Furthermore, the expression obtained 
for logistic grows permits discussion of the contribution of the shape of the limit cycle to 
stability. For t G [6i,T — 62], where 61^2 > are small, we find that 



< 7iW < 1 



(2^ 



which implies that the integral 7i(t)(it is small as well. Therefore limit cycles where the 
number of prey is extremely small for long times are not expected to be a basis for a stable 
ecosystem. For the AUee effect we find 



- 1] = ki{k2 - 7i)(7i - ^3) 
= A;i7i(/c2 + /cs - 271). 



1 + 71 



1 



1 



7i -h 7i + ^3 



(29) 



Considering the contribution of this term to the Floquet multiplier yields 
exp Ps[sx-l]dt^ = /i;i7i(fc2 + ^3 - 27i)(it^ 



exp ( I /ci7i(fc2 + k3)dt 



ki2-ffdt 



Increasing k2 and/or k^ will decrease stability. This is natural as these parameters represent 
the threshold to growth and the carrying capacity, providing another example for the paradox 
of enrichment. 

Note that the shape of the limit cycle can influence stability. In particular, the same 
conclusion to assumption (128|) holds. In the case of the AUee effect the de-stabilization effect 
for long periods of low prey density even enters quadratically in the term ki2'jfdt. This 
confirms the intuitive conclusion that imposing a threshold to growth is a de-stabilizing 
factor for non-equilibrium systems when the prey density is small. 

We proceed to consider the mortality of the predator. A very common functional form 
used in a large number of models is so-called density independent (linear) mortality 



M{Y) = kiY 



Pm\my - 1] = ki\my - 1] = A;i[l - 1] = 0. 
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Common Functional Forms 



G{X,Y) = kiXY (HoUing type I), 

G{X,Y) = 1^ (Hollingtypell), 

G{X,Y) = Ifg (Holling type III), 

G{X,Y) = Jl.^'Z.x^ (Holling type IV). 



Terms occurring in the Floquet multiplier 



(Holling 


type I) 


/3i[gx 


-1] = 








-1] = 


(Holling 


type II) 


Pi[gx 


-1] = 








-1] = 


(Holling 


type HI) 


I3i[9x 


-1] = 






(^2\gy 


-1] = 


(Holling 


type IV) 


HQx 


-1] = 






higy 


-1] = 



A;i72[l-1] = 0, 
A:i7i[l-1] = 0, 

fci72 r k2 _ 1] = fcl7l72 

fc2+7l '-fc2+7l J (fc2+7l)^' 

fcl7l72r_2fc2 _ 11 — fci7i72(fc2-7?) 

k2+li ^fc?+7? ^ (fc2+7?)2 • 



-1] 

[1 - 1] = 0, 



fcl7i 
fc2+7l 

A:i7i72 r 2A;2+7i _ 1 1 = 7172(^2-^37?) 
fc2+7i+fc37f '■fc2+7l+A:37i J (^2+71+^371 )^ ' 

fc2+7!+fc37t' ~ ^] = 0- 



Table 1: Calculation of the scale and elasticity functions for the predation term G{X^Y). 
The top panel lists four typical functional forms. In the bottom panel we calculate the terms 
occurring in the Floquet multiplier fl22|) . 



Therefore, linear predator mortality has no effect on the stability of the periodic solution. 

The interaction term between prey and predator is usually the most complicated and 
debated choice for the model. Some common choices are considered in Table [H The ob- 
servation that P2[gy — 1] vanishes for all functions considered in Table [T|, can be directly 
linked to the ecological assumption that predators hunt independently of each other. The 
functions that are therefore used in practice are generally linear in the density of predators 
and the impact of predator dependence on stability vanishes. The same assumption cannot 
generally be made for prey dependence of predation, leading to more complex expressions 
for the impact on stability. 

Therefore, we are going to make the assumptions 

= 1 and niy = 1 (30) 

from now. Regarding the Floquet multiplier formula ( |22l) the assumptions (|30ll simplify the 
situation to investigating 

A = exp(^^ /3,{s, - 1) - (3iig, - l)dty (31) 
The influence of and /3i on stability is not obvious since there is a non-trivial interaction 
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with the shape of the hmit cycle. The flow on moduh space given by (12^ simphfies to 

I3[ = /3i((/3.-/3i)^7.-(/3s-/3i) + (/52-/g™)), 

= ms- 1^1)9.. 

/?' = 

where we can view Pm as a parameter and simply drop the last equation. 




Figure 1: Dynamics in a specific example, (a) Stable periodic orbit 7(t) of (15^ (solid black) 
and two other trajectories (dashed magenta) with initial conditions marked by stars; the 
parameters are given in flMl) . Five points (black dots) are shown on the limit cycle for 
orientation purposes which are equally space over one period, (b) Scale functions in moduli 
space (black) for 7 solving fl32l) : a trajectory (solid magenta) with slightly perturbed initial 
conditions is also shown where the same elasticities as for the periodic orbit were used for 
numerical integration, (c) Time series of /3s (t) for part (b).(d) Time series 71 (solid black) 
and 72 (dashed black), (e) Scale functions associated to 7: Psif) (red), Pmit) (green), Piit) 
(blue) and /32(^) (cyan), (f) Elasticity function associated to 7: Sxif) (red), Qxit) (blue) and 
Qyit) (green); note that my = 1 = Qy. 



Example 5.1. For gaining an intuitive understanding one can consider the fiow on the 
moduli space in a specific example. The combination of logistic prey growth, Holling-type-II 
interaction and linear predator mortality gives us the Rosenzweig-MacArthur predator-prey 
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model that can produce periodic solutions 



(33) 



where we use the parameters 

fci = 2, k2 = 0.5, ks = 1, ki = 1, ^5 = 0.5. (34) 

Figure [1] shows that integrating a slightly perturbed initial condition trajectory does seem 
to diverge from the exact periodic solution in moduli space. Furthermore, even for a clas- 
sical planar predator-prey system, the scale and elasticity functions are quite complicated 
for non-equilibrium solutions. In fact, prescribing the elasticities is much more difficult than 
just picking a set of fixed parameters for equilibrium generalized models. 



To verify the necessary condition from Corollary 14.31 for periodic solutions we must ask 
for solvability of the boundary value problem (BVP) 

P[ = Pim-/3i)9.-{/3s-^i) + {^2-k,)), 

/3(0) = f3{T) forT>0, 

where P{t) := {I3s{t), Pi{t), P2{t))- It is well-known that BVPs can have one, many or no 
solutions [2]. Furthermore determining solvability conditions is usually not easy and even 
using numerical methods may be dangerous; for example, if a numerical algorithm fails to 
provide a solution to ( l35|) this may just be due to the numerical problems that can arise 
when solving BVPs [2]. 



6 Fourier Decomposition 

The previous discussion of specific functions and the Rosenzweig-MacArthur model motivate 
the need for a more concrete version of the moduli space conditions (1351) and of the Floquet 
multiplier f l3T|) . The natural step is to use a decomposition of the periodic functions into 
Fourier series; see Section 12.31 Using discrete convolution we can easily re- write the problem 
on moduli space. 

Proposition 6.1. Suppose we are given T-periodic elasticity functions s^, ruy, and gy. 
Then the Fourier coefficients of periodic scale functions have to satisfy the following set of 
algebraic equations 

^^^{k) = [/3,*(/3,-/3i)*(s,.-i)](fc), 

^-^L{k) = [/3„*(/32-/3™)*(^x.-i)](fc), , . 

^/3i(fc) = [/3i*((/3s-/3i)*^.. + (/32-/3™)*^7,-(4-/3i))](fc), ^ ' 

^-^m) = 02*{0s-Pl)*9. + 02-L)*9y-02-Pmmk), 

for all k & Z where we have also used the notation 1(0) = 1 and i(A;) = for k ^ and 
employed the obvious definition for addition of infinite sequences. 
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Proof. To complete the proof we only have to recall another basic fact from Fourier analysis. 
For two T-periodic functions we f , g we have 



oo oo 

2wi{k + m)t 
T 



fc=— oo m=—QO 

oo oo oo 

= E ^ hk)g{.n - k) = ^ {f * g){n) e" 



n=— oo fc=— oo 



This formula for Fourier coefficients of products of functions yields the right-hand side of 
equation f l36l) as a direct consequence of Theorem 14 .21 The left-hand side of equation f p6|) fol- 
lows from direct differentiation which is allowed since all our periodic functions are assumed 
to be sufficiently smooth; see Section 12.31 □ 
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Figure 2: Absolute value of the first nine Fourier coefficients < 4) associated to the 
stable periodic orbit •yit) of (15^ : the parameters are given in (IMI) . The coefficients of the 
phase space coordinates as well as the generalized elasticity and scale functions are shown. 



Since f l36|) is an infinite set of algebraic equations it may look like we have not consid- 
erably simplified the problem of finding scale functions that are consistent with prescribed 
elasticities. However, the rapid decay of Fourier coefficients provided by Theorem 12.11 al- 
lows us to approximate the solution of ( l36l) by focusing on the first few harmonics with 

\k\ < K, oo. 

Example 6.2 (Example 15.11 continued) . Just for illustration purposes we look a the Fourier 
coefficients of generalized scale and elasticity functions in an example. Figure [2] shows the 
results for the Rosenzweig-MacArthur model from Section |5] with k = 4. We can clearly see 
that the Fourier coefficients decay very rapidly; it is also interesting to observe that Sx:{k) is 
bimodal for logistic growth whereas the other coefficients show a uni-modal distribution for 
the first few harmonics. For the Rosenzweig-MacArthur model the algebraic relations (!36|) 
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on the Fourier coefficients become 



^/9i(A;) = [Pi*{0s-^)*9. + 02-L)-0s-mik), (37) 
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Figure 3: Absolute value of the first seventeen Fourier coefficients {\k\ < 17) for the left-hand 
and right-hand sides of the algebraic conditions (157|1 ; parameter values used are given in (IMj) . 
The black coefficients (lines shifted slightly left) are the coefficients of the derivatives (3'^, /3[ 
and P2 and the green coefficients (lines shifted slightly right) are associated to the periodic 
functions on the right-hand side of fl37j) . The agreement of the two sets of coefficients is 
clearly visible. 

Figure |3] shows the values of the Fourier coefficients for the Rosenzweig-MacArthur ex- 
ample where we see that the algebraic conditions f p7|) are satisfied as proven in Proposition 
16.11 Furthermore, it is evident that due to the convolution a wider support km is necessary 
i.e. the algebraic equations (136!) must be satisfied for \k\ < km where km > 1^ and k is our 
truncation for the Fourier coefficients of the phase space periodic orbit. 



7 Stability Analysis 

To get a better understanding of stability we can also use the Fourier series approach to 
re-express the Floquet multiplier (|3T1) given by 

A = exp (^j^ /3,(s, - 1) - Pi{g, - 1) + P2{gy - 1) - Pm{my - l)dt^ . (38) 
The next results shows how the different Fourier coefficients enter in formula fl38l). 



Theorem 7.1. For the non- equilibrium generalized predator-prey model with Qy = 1 = ruy 
the single Floquet multiplier of a T -periodic orbit is given by 

A = exp I T I [ 4*(g..^- i)m - 0i*{g^- 1)]( 0) | | = exp(T(Ci - C^)) (39) 

= :Ci =:C2 
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i.e. whether |A| > 1 or |A| < 1 depends only on the difference of two zeroth-order Fourier 
coefficients Ci and C2 that arise from two discrete convolutions. 

Proof. We start by looking at the first summand in integral in (138|1 which gives 
/ Psis. - l)dt = / V * (s, - i)]{k) e^^'^'/^dt 

k=-oo 

00 „T 

= Yl 0s*is.-l)m / e'^'^'^'/^dt 

k=-oo 

_ r for ^ 

~ \ t [4*(s,-i)](0) for fc = 

where the last step follows from the fact that e'^'^^^^^^dt = for k ^ 0. From this 
calculation we find Ci and in a similar way also C2. Using the two factors Ci^2 and rriy = 
1 = iia fl38l) the result fl39|) follows. Since T > the modulus |A| only depends on the 
difference of Ci and C2; if Ci — C2 > then |A| > 1 and if Ci — C2 < we obtain |A| < 1. □ 

Before we consider in more detail the dependency of stability on Ci and C2 we briefly 
investigate the influence of the period T. Although T does not effect the stability of a periodic 
orbit directly it does have an interesting biological interpretation. If T ^ 1 then the period 
amplifies stability and instability. For example, when Ci — C2 > then a long period moves 
the multiplier even further away from |A| = 1 and trajectories near the unstable periodic 
orbit will escape vary quickly. On the other hand, if Ci — C2 < and |A| < 1 then a very 
large period T moves the multiplier even closer to the super- attracting case < |A| ^ 1. A 
very short period < T ^ 1 has the effect of moving the multiplier very close to |A| ~ 1. 
This means that when the periodic orbit is unstable, it will take a very long time to escape 
from it. The last effect can be interpreted as inducing met a- stability i.e. when the period of 
the predator-prey cycle is short then the predator-prey system stays near a metastable state 
for a long time although it is eventually unstable. This could lead to the conjecture that 
fast oscillations could be beneficial to survival for predator-prey populations during periods 
when external parameters entering Ci and C2 drive the system, potentially only temporarily, 
to a state when |A| > 1. 

As a next step, we want to understand better how the Fourier coefficients of /3s, 
and influence Ci and C2. 

Proposition 7.2. The two constant Ci^2 o^^e given by 

00 

Ci = ^s{0){U0)-l) + 2j2iM^ik)]Re[s,{k)] + Im[^smM^xm) (40) 

k=l 

00 

C2 = PmiU0)-l) + 2j2iMPi{k)]Re[Uk)] + Irnmk)^^^^ (41) 

k=l 
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Proof. Given two sequences {/(fc)}fcL-oo {9{k)}T=-oo '^^ Fourier coefficients for two real- 
valued functions a direct calculation yields 

oo 

(/*^)(o) = hm-k) 

= f{o)m + J2 fx-mk) + E hk)gi-k) 

A;>0 fc>0 

where we have used f{—k) = f{k) and the real-valuedness f = f, g = g in the last step. 
Next, observe that 

Wmk) + f{k)W) = 2{Re[f{k)]Re[g{k)] + Im[/>)]Im[^(A;)]). 

Now the formulas (140|) - (14T|) follow immediately as f3s(t), f3i(t), Sx{t) and gxit) are all real- 
valued. □ 

In practice, we never use the infinite sum formulas from Proposition 17.21 but truncate 
them at a finite order. Using the explicit formulas for Ci and C2 we can directly draw 
several conclusions regarding periodic solutions depending on generalized scale and elasticity 
functions (recall: we still use gy = 1 = my). If all Fourier coefficients of higher-order k >1 
are small, then stability of periodic solutions is dominated by the terms 

Ci ^ 4(0)(S,.(0) - 1) and ~ PiimM - !)• 

Since the scale functions are always positive the time averages of the elasticity functions 
Sa;(0) and cjxiS^) determine the signs of Ci and C2. Therefore, average sub-linear elasticity 
"Sa;(0) < 1 and average super-linear conversion ^^^(O) > 1 enhance stability. In ecological 
terms < s^(0) < 1 means that, on average, the prey growth should be limited by external 
factors and S2:(0) > 1 means that, on average, the predation rate should be sensitive to 
prey abundance; see also [21] for an interpretation of the generalized parameters for the 
equilibrium case. Both conditions make intuitive sense: if the prey grows without external 
limitation then solutions may be expected to diverge from a periodic solution and become 
unbounded while insensitivity of predation to prey growth could potentially drive a system to 
extinction. Of course, also the inverse relationships hold so that Sa:(0) > 1 and ^x(O) < 1 act 
towards de-stabilization. In this context, the scale functions act as amplifiers. For example, 
if Ci < and C2 > then a large average growth rate /3s (0) ^ 1 and a large conversion rate 
/3i(0) ^ 1 will enhance stability even more since the Floquet multiplier moves closer to the 
super- attracting regime A ~ 0. In this case, initial conditions will be attracted much quicker 
to a stable limit cycle. If 

4(0)(s,.(0) - 1) - /3i (0)^(0) -1)^0 i.e. 4(0)(s,.(0) - 1) ^ /3i (0)^(0) - 1) 

the leading-order terms between growth and predation balance and the stability properties 
are dominated by higher-order harmonics. The leading-order terms also become irrelevant 
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(a) (b) (c) 




-2-1012 -2-1012 -2-1012 



Re Re Re 

Figure 4: Illustration how the location of Fourier coefficients for generalized elasticity and 
scale functions influence stability. Here we focus on the first two higher-order harmonics 
(ffist coefficient = solid line, second coefficient = dashed line) of f3s(t) (red) and Sx(t) (blue) 
which influence the term Ci in Theorem 17. 1[ (a) vr-phase shift gives Ci < 0. (b) Small phase 
shift gives Ci > 0. (c) Competition between ffist- and second-order harmonics. 

for elasticity functions which average close to one 

S,(0) ^ 1 and ^,(0) ^ 1. 

In this scenario we have to focus on the relations between the higher-order Fourier coeffi- 
cients of Ps and as well as Pi and g^. Let us assume for simplicity that 5^.(0) = 1 = ^'^(0) 
so that we can focus on the higher harmonics. Then stability enhancing conditions are 

oo 

Ci = 5^(Re[/3,(fc)]Re[S,(fc)] + Im[4(A;)]Im[S,(fc)]) <0, 

k=l 
oo 

C2 = ^(Re[/3i(fc)]Re[^,(fc)]+Im[/3i(fc)]Im[^,(A:)]) >0. 

k=l 

Figure m depicts several different situations in the complex plane for the ffist two higher-order 
harmonics of Ps{k) and Sa;{k) {k = 1,2). In Figure |l](a) the ffist two higher-harmonics are 
in "anti-phase" so that the angles between the coefficients are separated by vr. This means 
that 

Re[Ps{k)]Re[s^{k)] < and Im[/3,(fc)]Im[s^(fc)] < 

for = 1, 2. In such a situation, we expect that Ci < by disregarding higher orders so 
that stability is enhanced. 

Figure IH^b) shows the situation where there is only a small phase difference between the 
coefficients ("in-phase") which gives 

Re[Ps{k)]Re[s^{k)] > and lm[Psik)]lm[s^{k)] > 0. 

There is also a possible situation where a competition between the different order harmonics 
arises as illustrated in Figure |l](c). We can now also give an ecological interpretation of these 



21 



conditions. Stability Ci < is enhanced if Sx{t) and (3s{t) oscillate with a phase separation 
near tt which means that a period of high sensitivity of prey abundance should coincide with 
a period of low prey growth and vice versa. Note that these conditions also make sense 
intuitively and suggest that prey growth is most efficient if there is a small number of prey 
and there are no limiting factors from the environment. Similar considerations also apply to 
the stability enhancing condition C2 > 0. A small phase separation between Pi{t) and gx{t) 
increases stability of the predator-prey limit cycle. Observe that gx{t) can be interpreted as 
the dependence of predation on prey abundance and Pi{t) as a predation rate (normalized 
by the total number of prey) |2l]; the stability conditions mean that a high predation rate 
should coincide with a high dependence of predation on prey abundance. In other words, 
if the dominating factor to prey abundance is predation then it is good for the predator to 
hunt a lot to increase stability of the limit cycle. 

Note that although the conclusions stated above seem to be "obvious" in an ecological 
context, it is by no means clear how to prove them. That they can be obtained by an analysis 
of nonlocal generalized models underlines the applicability of the approach. 

8 Sampling 

Recall that due to Proposition 14.11 it was straightforward for equilibrium generalized models 
to choose a set of generalized parameters, just random sampling produces a set of parameters 
that is consistent with at least one specific model. Random sampling of generalized param- 
eters has been exploited to correlate different aspects of the dynamical system to stability 
[ISlllG]. For non-equilibrium systems we must certainly check the necessary condition from 
Corollary 14. 3[ One possibility is the following algorithm which allows sampling of elasticity 
and scale functions: 

(Al) Prescribe a set of T-periodic elasticity functions by their Fourier coefficients. For 
simplicity we will always choose T = 1 and assume Qy = 1 = rriy. 

(A2) Choose a truncation order for the algebraic system (1371) so that the necessary 
condition reads 

= -2mk^sik) + 0s*0s-^i)*is.-i)]ik) =: c,(A;), 

= -27izkP2{k) + [P2*m-^)*g.)m =■■ C2{k), 

for \k\ < Km- 

(A3) Define a new variable that collects all the Fourier coefficient values for (3s, Pi, P2, Pm, 
Sx and gx 

X := (/3,(0), /3,(1), . . . , PsM, /3i(0), . . . , . . .) e C^^^-^+i) ^ ^i^i-M+i) 

where x contains all the information about the scale and elasticity functions since the 
negative index coefficients can be obtained by complex conjugation. 
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(A4) Define a function 

F{X) := ||Re(c.)f + ||Re(cOf + ||Re(c2)f + ||Im(c,)f + HMcOf + WMcW 

where we view c^, ci and C2 as vectors of dimension 2k,m + 1 and real and imaginary 
parts are applied component-wise. 

(A5) Observe that F{Xq) = if and only if the Fourier coefficients encoded in Xq satisfy 
the algebraic equations in (A2). Therefore we can attempt to solve the following 
optimization problem 

Xm ■■= min{F(X) : X G m12(«m+i)| (42) 
with a random initial condition, say x = xi. 

Solving the optimization problem for different random initial conditions is expected to 
yield different values for Xm that solve the algebraic constraint in (A2). This means that we 
get a set of Fourier coefficients {Xm{l)}f=i where L denotes the sample size and the index 
/ G N indicates the dependence on the initial condition. 

The main technical difficulty of the algorithm (Al)-(A5) is that it involves the solution 
of the optimization problem fH2]) . This is computationally much more expensive than the 
direct random sampling for equilibrium generalized models. It is known [H] that the main 
computational cost in optimization is often given by the difficulty of the function evaluations 
of F{x). For our case, this seems to be the case since we have to compute several discrete 
convolutions to evaluate F{x). However, the convolution computation is inexpensive due to 
the Fast Fourier Transform [31j. 

Now we want to demonstrate that the algorithm can be used for a sampling analysis of 
stability similar to the one used in [22]. Let us point out that we do not attempt a full 
detailed statistical analysis here but that we only aim at a proof-of-principle. We solved 
fH2]) for 110000 initial conditions for km = 2 using a standard algorithm for nonlinear 
optimization [521 ES] • Each sequence of Fourier coefficients in the initial condition consists 
of five real numbers e.g. 

/3,(0), Re(/3,(1)), Im(4(l)), Re(/3,(2)), Im(/3,(2)), (43) 

which were sampled uniformly and independently from the interval [0.5, 1.5]. We discarded 
all solutions of the optimization algorithm that did not satisfy the positivity condition 

4(0) > 0, /3i(0) > 0, /32(0) > 0, /3„(0) > 0. 

which is required by the definition of the scale functions and the invariance of the positive 
quadrant for the moduli space flow. The 63587 remaining solutions Xm(/) satisfled the the 
optimization problem (and therefore the moduli space flow) at least up to a tolerance of 
10~^ i.e. |a;m(OI < for all /; the average value was E[a;m(0] ~ " 10~^. We have also 
calculated the single Floquet multiplier A; associated to each solution using Proposition 17.21 
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Figure 5: Histogram of the 5-63587 Fourier coefficients for fSg obtained from the optimization 
of (H2|) with uniformly sampled initial conditions ( 143|) . The columns show the five different 
real numbers with their observed number on the vertical axes. The ffist row shows coefficients 
associated to a stable Floquet multiplier and the second row those with an unstable Floquet 
multiplier. Observe that the number of stable coefficients is substantially larger than the 
number of unstable ones. 



Figure [5] shows some of the output of the computation. We plot the Fourier coefficients 
associated to the scale function Ps- The top row in Figure |5] corresponds to coefficients with 
stable periodic orbit (|A| < 1) and the bottom row to coefficients with an unstable periodic 
orbit (|A| > 1). We see that, despite the initial uniform sampling, the results for each 
coefficient of (3s closely resemble normal distributions. The same observation also applies for 
the other scale and elasticity functions. In total we find that 37873 solutions associated to a 
stable multiplier and 25714 unstable ones. From this discrepancy one may either conjecture 
that the moduli space flow constraint could bias ecosystem towards stabihty or that our 
choice of initial uniform random sampling over a particular region in parameter space causes 
the bias towards stability. 

In Table |2] we list mean and variance for each coefficient. Several observations can be 
made based on Table El The scale functions (32 and have a much bigger variance than 
(3s and (3i. This could indicate that the prey growth rate and the prey-per-capita predation 
rate have to obey much smaller ranges in ecosystems compared to the predator-per-capita 
rates describing consumption and mortality. It is also interesting that the mortality rate 
(3m allows for much larger amplitude higher-order harmonics whereas e.g. |/3s(2)| is always 
comparatively small. The elasticities show no consistent variance decay towards higher- 
order harmonics although the coefficients themselves seem to decay. From the ecological 
perspective this suggest that predator-prey systems may exhibit a wide diversity in terms of 
sensitivities and g^. 
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/3.(0) 


Rc(/§41)) 


Im(/3s(l)) 


Re(/§42)) 


Im(/§42)) 


mean (stable) 
variance (stable) 
mean (unstable) 
variance (unstable) 


1.6021 
0.8840 
1.4409 
0.9129 


0.2207 
0.3226 
0.1936 
0.2863 


0.5035 
0.4000 
0.3863 
0.3409 


-0.0117 
0.1000 

-0.0232 
0.1014 


-0.1205 
0.1202 

-0.0591 
0.1381 




^i(O) 


Re(/3i(l)) 


M^i(i)) 


Re(^i(2)) 


Im(/3i(2)) 


mean (stable) 
variance (stable) 
mean (unstable) 
variance (unstable) 


1.3182 
0.6575 
1.4502 
0.7578 


0.3275 
0.3050 
0.1778 
0.3213 


0.3219 
0.2431 
0.3445 
0.2796 


0.0157 
0.1502 
-0.0686 
0.1679 


0.2875 
0.1226 
0.2649 
0.1344 




/32(0) 


Re(/32(1)) 


Im(fe(l)) 


Re(fe(2)) 


Im(/32(2)) 


mean (stable) 
variance (stable) 
mean (unstable) 
variance (unstable) 


2.0805 
1.8865 
1.8914 
1.8228 


0.4823 
0.6716 
0.3399 
0.5437 


0.5761 
0.7486 
0.3835 
0.5802 


-0.0583 
0.2135 

-0.0375 
0.1593 


0.2154 
0.2970 
0.0502 
0.2070 


/3m (fc) 


An(0) 


Re(/§m(l)) 


Im(/3m(l)) 


Re(^m(2)) 


Im(/3m(2)) 


mean (stable) 
variance (stable) 
mean (unstable) 
variance (unstable) 


1.6184 
1.3568 
1.7736 
1.5189 


0.9865 
22090 
1.0465 
2.2642 


0.4756 
2.3338 
0.6171 
2.7503 


1.7548 
2.9748 
1.7785 
2.9247 


0.4965 
3.1220 
0.7635 
3.4595 




s^(0) 


Re(s,(l)) 


Im(s,(l)) 


Re(s,(2)) 


Im(s^(2)) 


mean (stable) 
variance (stable) 
mean (unstable) 
variance (unstable) 


1.5988 
2.2079 
2.5967 
3.7412 


1.0598 
2.8519 
1.5343 
3.0413 


1.6099 
2.2995 
1.9850 
3.5637 


1.5222 
2.4361 
1.5559 
2.7521 


0.8393 
3.1797 
1.2697 
2.9364 


9x{k) 


MO) 


Rc(§4i)) 


Im(34l)) 


Re(s42)) 


Im(g,(2)) 


mean (stable) 
variance (stable) 
mean (unstable) 
variance (unstable) 


2.7354 
4.3300 
1.4787 
2.3056 


1.9722 
3.6009 
1.6774 
3.4147 


2.4554 
3.4094 
1.8789 
3.4502 


0.9165 
2.9125 
1.2302 
3.5286 


1.3490 
2.6612 
1.1109 
3.1522 



Table 2: Mean and variance for the Fourier coefficients obtained from optimization (solution 
of the moduli space flow). Coefficients associated to stable and unstable Floquet multipliers 
are considered separately. 



To understand how the different coefficients relate to stability we calculate the Pearson 
correlation coefficient. For two vectors of observations {a/} and {bi} it is defined as 

. ^^ Eiiai-E[a])ik-E[b]) 

Figure |6] shows r(a, A;) where a is a sequence of real or imaginary parts of the Fourier 
coefficients e.g. {ai} = {Iie{P s,i{k))}- One important conclusion to draw from the correlation 
coefficients is that although a Fourier coefficient does not appear in the stability formula 
for the Floquet multiplier it may still correlate positively or negatively with stability. For 
example, /32(1) and /32(2) show a negative correlation with Floquet multiplier. This effect 
can be caused by the fact that the scale and elasticity functions are not independent i.e. 
they are related via the moduli space flow. 

It is very important to observe that we can recover conclusions, which we found already 
analytically in Section [TJ from the statistical analysis. For example, the coefficient ^^^(O) 
correlates negatively with stability which means that decreasing it increases the Floquet 
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Figure 6: Pearson correlation coefficient r = r{a,b) between stability and the different 
harmonics; positive correlation is indicated in blue and negative correlation in black. Each 
panel represents correlation for the Fourier coefficients from left to right. For example, the top 
left panel shows the values r{Xi, 4/0)), r{Xi, Re0s,ii^))), r(A,, Im(/3,/l))), r{Xi, Re(4/2))), 
r(A;, Im(/^s^;(2))) from left to right where A; is the Floquet multiplier with index /. 



multiplier and acts towards destabilization. This is precisely the result we have already 
obtained analytically in Section [3 Let us point out again that the basic statistical analysis 
we have provided is incomplete but that it definitely does show that the proposed sampling 
techniques based on the FFT, optimization and correlations can help to understand stability 
of periodic solutions. 



9 Outlook 

In this paper we have extended the method of generalized modeling from equilibrium to 
non-equilibrium systems. This extension has been achieved in the context of a classical 
predator-prey system with periodic solutions. The main re-scalings and definitions from the 
equilibrium case can be carried over to periodic orbits. However, the resulting generalized 
ODEs differ from the steady state case in several respects. The algebraic form is different 
due to the time dependent re-scalings and also the generalized parameters become time- 
dependent elasticity and scale functions. The Jacobian A{t) of the system has to be analyzed 
using Floquet theory that describes the stability of periodic orbits. For planar vector fields 
we have been able to use Liouville's formula 

A = exp (^j^ Tr{A{t))d?j 

which facilitated several analytical calculations. We have discovered that the generalized 
elasticity and scale functions have to satisfy a flow moduli space. Then we used Fourier 
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analysis to find computable conditions from the moduli flow. Discrete convolutions turned 
out to be the key to stability analysis providing explicit interpretable stability results. In the 
last part of the paper, we suggested a sampling algorithm that uses optimization methods to 
find elasticity and scale functions that satisfy the (algebraic) moduli space flow. During our 
analysis we have also obtained several ecological conclusions about arbitrary predator-prey 
models that can be written in the generalized form ([3]). 

In principle, we can extend the theory described here without any technical problems 
to limit cycles in A^-dimensional systems for N > 2; see also ^33j regarding generalizations 
to in the equilibrium context. The main difference in will be that we have to 
compute several the Floquet multipliers numerically since Liouville's formula only provides 
the product of the eigenvalues. 

One can also consider a generalization to non-equilibrium system beyond periodic orbits. 
For example, the generalized model fl21l) . as well as Theorem 14.21 on the moduli space flow, 
carry over directly to other situations such as homoclinic trajectories [33] or chaotic dynamics 
[23] . For instance, instead of posing periodic boundary conditions of the form 



we have to impose other conditions on the solution of the moduli space flow. For homoclinic 
orbits we need the boundary conditions 



i.e. that we have asymptotic limits of the generalized elasticity and scale functions to their 
value at a saddle-point equilibrium. For chaotic dynamics one must search for aperiodic 
bounded trajectories in moduli space. Note that this raises interesting mathematical as well 
as application questions. For example, the moduli space flow may provide new insights when 
a dynamical system may be chaotic. Finally, also our sampling analysis can obviously be ex- 
tended. Beyond a more detailed statistical validation, we could consider higher-dimensional 
food webs [22] which leads to a problem in R^. 
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